      SUBROUTINE STEP_NCYC (ISTEP,DT,ALPH)
C--
C--   SUBROUTINE STEP_NCYC (ISTEP,DT,ALPH)
C--
C--   author : Daisuke Hotta (DH)
C--   Purpose: perform one time step by Lorenz 4-cycle
C--            using the following scheme:
C--   w    = wlist(mod(ISTEP,Nperiod))
C--   F(2) = w * [ T_dyn(F(1)) + T_phy(F(1)) ] + (1-w) * F(2)
C--   F(1) = F(1) + DT*F(2)
C--   where F(2) is the accumulated tendency and
C--   wlist is an array of weights defined in "ncycparm.h"
C--   Input: 
C--   ISTEP      : time step index (This determines the weight)
C--   DT = time step (if DT < or = 0, tendencies are computed but 
C--                   no time stepping is performed)
C--   ALPH = 0   : forward step for gravity wave terms
C--   ALPH = 1   : backward implicit step for g.w.
C--   ALPH = 0.5 : centered implicit step for g.w.
C-- 
C--   Modified common blocks : DYNSP1, DYNSP2
C--
      include "atparam.h"
      include "atparam1.h"

      include "com_dyncon0.h"
      include "com_hdifcon.h"

      include "com_dynvar.h"

      COMPLEX VORDT(MX,NX,KX), DIVDT(MX,NX,KX), TDT(MX,NX,KX),
     *        PSDT(MX,NX), TRDT(MX,NX,KX,NTR)

      COMPLEX CTMP(MX,NX,KX)
      COMPLEX ZERO
      include "ncycparm.h"

      ZERO = (0.,0.)

      iitest=1
      if(iitest.eq.1) print*, ' inside step'

C--   1. Computation of grid-point tendencies
C         (converted to spectral at the end of GRTEND)

      if (iitest.eq.1) print*,' call grtend'
      CALL GRTEND (VORDT,DIVDT,TDT,PSDT,TRDT,1,1)

C--   2. In this scheme, spectral tendencies are not computed every time step
        
C--   2.5 Update accumulated tendencies and overwrite tendencies with
C--       the accumulated tendencies
      CALL UPD_ACM_TND(ISTEP, 1, PS, PSDT)
      
      CALL UPD_ACM_TND(ISTEP,KX,VOR,VORDT)
      CALL UPD_ACM_TND(ISTEP,KX,DIV,DIVDT)
      CALL UPD_ACM_TND(ISTEP,KX,  T,  TDT)
      
      DO ITR=1,NTR
         CALL UPD_ACM_TND(ISTEP,KX,
     &        TR(1,1,1,1,ITR),  TRDT(1,1,1,ITR))
      ENDDO
      
C--   3. Horizontal diffusion
      
      if (iitest.eq.1) print*, ' biharmonic damping '
      
C     3.1 Diffusion of wind and temperature
      
      CALL HORDIF (KX,VOR,VORDT,DMP, DMP1)
      CALL HORDIF (KX,DIV,DIVDT,DMPD,DMP1D)
      
      DO K=1,KX
         DO M=1,MXNX
            CTMP(M,1,K)=T(M,1,K,1)+TCORH(M,1)*TCORV(K)
         ENDDO
      ENDDO
      
      CALL HORDIF (KX,CTMP,TDT,DMP,DMP1)
      
C     3.2 Stratospheric diffusion and zonal wind damping
      
      SDRAG=1./(TDRS*3600.)
      DO N=1,NX
         VORDT(1,N,1)=VORDT(1,N,1)-SDRAG*VOR(1,N,1,1)
         DIVDT(1,N,1)=DIVDT(1,N,1)-SDRAG*DIV(1,N,1,1)
      ENDDO
      
      CALL HORDIF (1,VOR, VORDT,DMPS,DMP1S)
      CALL HORDIF (1,DIV, DIVDT,DMPS,DMP1S)
      CALL HORDIF (1,CTMP,TDT,  DMPS,DMP1S)
      
C     3.3 Chech for eddy kinetic energy growth rate 
      
C     CALL CGRATE (VOR,DIV,VORDT,DIVDT)
      
C     3.4 Diffusion of tracers
      
      DO K=1,KX
         DO M=1,MXNX
            CTMP(M,1,K)=TR(M,1,K,1,1)+QCORH(M,1)*QCORV(K)
         ENDDO
      ENDDO
      
      CALL HORDIF (KX,CTMP,TRDT,DMPD,DMP1D)
      
      IF (NTR.GT.1) THEN
         DO ITR=2,NTR
            CALL HORDIF (KX,TR(1,1,1,1,ITR),TRDT(1,1,1,ITR),
     &           DMP,DMP1)
         ENDDO
      ENDIF
      
C--   4. Time integration using Lorenz N-cycle

      IF (DT.LE.0.) RETURN
      
      if (iitest.eq.1) print*,' time integration'
      
      CALL TIMINT (DT,1,PS,PSDT)
      
      CALL TIMINT (DT,KX,VOR,VORDT)
      CALL TIMINT (DT,KX,DIV,DIVDT)
      CALL TIMINT (DT,KX,T,  TDT)

      DO ITR=1,NTR
        CALL TIMINT (DT,KX,TR(1,1,1,1,ITR),TRDT(1,1,1,ITR))
      ENDDO
C-- 5. Apply Semi-implicit correction 
      IF ( MOD(ISTEP,Ncyc) .eq. (Ncyc-1) ) THEN
c     clear tendencies
         DO K=1,KX
            DO M=1,MXNX
               DIVDT(M,1,K)=ZERO
               TDT(M,1,K)  =ZERO
            ENDDO
         ENDDO
         DO M=1,MXNX
            PSDT(M,1) = ZERO
         END DO
c     compute spectral tendencies
      if (iitest.eq.1) print*,' call sptend'
      CALL SPTEND (DIVDT,TDT,PSDT,1)

c     implicit correction 
         if (iitest.eq.1) print*,' call implic'
         CALL IMPLIC (DIVDT,TDT,PSDT)
         print *, MAXVAL(abs(DIVDT))
         print *, MAXVAL(abs(TDT))
         print *, MAXVAL(abs(PSDT))
c     apply the correction
         CALL TIMINT (Ncyc*DT,1,PS,PSDT)
         CALL TIMINT (Ncyc*DT,KX,DIV,DIVDT)
         CALL TIMINT (Ncyc*DT,KX,T,  TDT)
         END IF

      RETURN
      END   

      SUBROUTINE HORDIF (NLEV,FIELD,FDT,DMP,DMP1)
C--
C--   Aux. subr. HORDIF (NLEV,FIELD,FDT,DMP,DMP1)
C--   Purpose : Add horizontal diffusion tendency of FIELD 
C--             to spectral tendency FDT at NLEV levels
C--             using damping coefficients DMP and DMP1
C--
      include "atparam.h"

      COMPLEX FIELD(MXNX,NLEV), FDT(MXNX,NLEV)
      REAL    DMP(MXNX), DMP1(MXNX)

      DO K=1,NLEV
        DO M=1,MXNX
          FDT(M,K)=(FDT(M,K)-DMP(M)*FIELD(M,K))*DMP1(M)
        ENDDO
      ENDDO

      RETURN
      END

      SUBROUTINE TIMINT (DT,NLEV,FIELD,FDT)
C--
C--   Aux. subr. TIMINT (DT,NLEV,FIELD,FDT)
C--   Purpose : Perform time integration of FIELD at NLEV levels
C--             using tendency FDT
C--
      include "atparam.h"

      COMPLEX FIELD(MXNX,NLEV,2), FDT(MXNX,NLEV), FNEW(MXNX)

      IF (IX.EQ.IY*4) THEN
        DO K=1,NLEV
          CALL TRUNCT (FDT(1,K))
        ENDDO
      ENDIF
      
      DO K=1,NLEV
         DO M=1,MXNX
            FIELD(M,K,1) = FIELD(M,K,1) + DT*FDT(M,K)
         ENDDO
      ENDDO
      
      RETURN
      END

      SUBROUTINE UPD_ACM_TND (ISTEP,NLEV,FIELD,FDT)
C--
C--   Aux. subr. UPD_ACM_TND (ISTEP,NLEV,FIELD,FDT)
C--   Purpose : Update accumulated tendency
C--             
C--
      include "atparam.h"
      INTEGER ISTEP
      COMPLEX FIELD(MXNX,NLEV,2), FDT(MXNX,NLEV)
      REAL W
      include "ncycparm.h"

      W = WLIST( MOD(ISTEP, Nperiod)+1 )
      DO K=1,NLEV
         DO M=1,MXNX
            FIELD(M,K,2) = W*FDT(M,K) + (1.-W)*FIELD(M,K,2)
            FDT(M,K) = FIELD(M,K,2)
         ENDDO
      ENDDO
      print *, 'DBG: in UPD_ACM_TND: maxval(FDT)=', 
     &  maxval(abs(FDT))
      RETURN
      END

      SUBROUTINE CGRATE (VOR,DIV,VORDT,DIVDT)
C--
C--   SUBROUTINE CGRATE (VOR,DIV,VORDT,DIVDT)
C--
C--   Purpose: Check growth rate of eddy kin. energy 
C--   Input  : VOR    = vorticity
C--            DIV    = divergence
C--            VORDT  = time derivative of VOR
C--            DIVDT  = time derivative of DIV
C--
      include "atparam.h"
      include "atparam1.h"

      COMPLEX VOR(MX,NX,KX), VORDT(MX,NX,KX), 
     &        DIV(MX,NX,KX), DIVDT(MX,NX,KX), TEMP(MX,NX)

      GRMAX=0.2/(86400.*2.)

      CDAMP=0.

      DO K=2,KX

        GRATE=0.
        RNORM=0.

        CALL INVLAP (VOR(1,1,K),TEMP)

        DO N=1,NX
         DO M=2,MX
           GRATE=GRATE-REAL(VORDT(M,N,K)*CONJG(TEMP(M,N)))
           RNORM=RNORM-REAL(  VOR(M,N,K)*CONJG(TEMP(M,N)))
         ENDDO
        ENDDO

        IF (GRATE.GT.GRMAX*RNORM) CDAMP =
     &      MAX(CDAMP,0.8*GRATE/RNORM)
C    &      MAX(CDAMP,(GRATE*GRATE)/(GRMAX*RNORM*RNORM))

      ENDDO

      IF (CDAMP.GT.0.) THEN

        print *, ' rot. wind damping enabled'

        DO K=1,KX
          DO N=1,NX
           DO M=2,MX
             VORDT(M,N,K)=VORDT(M,N,K)-CDAMP*VOR(M,N,K)
           ENDDO
          ENDDO
        ENDDO

      ENDIF


      CDAMP=0.

      DO K=2,KX

        GRATE=0.
        RNORM=0.

        CALL INVLAP (DIV(1,1,K),TEMP)

        DO N=1,NX
         DO M=2,MX
           GRATE=GRATE-REAL(DIVDT(M,N,K)*CONJG(TEMP(M,N)))
           RNORM=RNORM-REAL(  DIV(M,N,K)*CONJG(TEMP(M,N)))
         ENDDO
        ENDDO

        IF (GRATE.GT.GRMAX*RNORM) CDAMP =
     &      MAX(CDAMP,0.8*GRATE/RNORM)
C    &      MAX(CDAMP,(GRATE*GRATE)/(GRMAX*RNORM*RNORM))

      ENDDO

      IF (CDAMP.GT.0.) THEN

        print *, ' div. wind damping enabled'

        DO K=1,KX
          DO N=1,NX
           DO M=2,MX
             DIVDT(M,N,K)=DIVDT(M,N,K)-CDAMP*DIV(M,N,K)
           ENDDO
          ENDDO
        ENDDO

      ENDIF

C--
      RETURN
      END
